{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "<style type='text/css'>\n",
       ".hll { background-color: #ffffcc }\n",
       ".c { color: #408080; font-style: italic } /* Comment */\n",
       ".err { border: 1px solid #FF0000 } /* Error */\n",
       ".k { color: #008000; font-weight: bold } /* Keyword */\n",
       ".o { color: #666666 } /* Operator */\n",
       ".ch { color: #408080; font-style: italic } /* Comment.Hashbang */\n",
       ".cm { color: #408080; font-style: italic } /* Comment.Multiline */\n",
       ".cp { color: #BC7A00 } /* Comment.Preproc */\n",
       ".cpf { color: #408080; font-style: italic } /* Comment.PreprocFile */\n",
       ".c1 { color: #408080; font-style: italic } /* Comment.Single */\n",
       ".cs { color: #408080; font-style: italic } /* Comment.Special */\n",
       ".gd { color: #A00000 } /* Generic.Deleted */\n",
       ".ge { font-style: italic } /* Generic.Emph */\n",
       ".gr { color: #FF0000 } /* Generic.Error */\n",
       ".gh { color: #000080; font-weight: bold } /* Generic.Heading */\n",
       ".gi { color: #00A000 } /* Generic.Inserted */\n",
       ".go { color: #888888 } /* Generic.Output */\n",
       ".gp { color: #000080; font-weight: bold } /* Generic.Prompt */\n",
       ".gs { font-weight: bold } /* Generic.Strong */\n",
       ".gu { color: #800080; font-weight: bold } /* Generic.Subheading */\n",
       ".gt { color: #0044DD } /* Generic.Traceback */\n",
       ".kc { color: #008000; font-weight: bold } /* Keyword.Constant */\n",
       ".kd { color: #008000; font-weight: bold } /* Keyword.Declaration */\n",
       ".kn { color: #008000; font-weight: bold } /* Keyword.Namespace */\n",
       ".kp { color: #008000 } /* Keyword.Pseudo */\n",
       ".kr { color: #008000; font-weight: bold } /* Keyword.Reserved */\n",
       ".kt { color: #B00040 } /* Keyword.Type */\n",
       ".m { color: #666666 } /* Literal.Number */\n",
       ".s { color: #BA2121 } /* Literal.String */\n",
       ".na { color: #7D9029 } /* Name.Attribute */\n",
       ".nb { color: #008000 } /* Name.Builtin */\n",
       ".nc { color: #0000FF; font-weight: bold } /* Name.Class */\n",
       ".no { color: #880000 } /* Name.Constant */\n",
       ".nd { color: #AA22FF } /* Name.Decorator */\n",
       ".ni { color: #999999; font-weight: bold } /* Name.Entity */\n",
       ".ne { color: #D2413A; font-weight: bold } /* Name.Exception */\n",
       ".nf { color: #0000FF } /* Name.Function */\n",
       ".nl { color: #A0A000 } /* Name.Label */\n",
       ".nn { color: #0000FF; font-weight: bold } /* Name.Namespace */\n",
       ".nt { color: #008000; font-weight: bold } /* Name.Tag */\n",
       ".nv { color: #19177C } /* Name.Variable */\n",
       ".ow { color: #AA22FF; font-weight: bold } /* Operator.Word */\n",
       ".w { color: #bbbbbb } /* Text.Whitespace */\n",
       ".mb { color: #666666 } /* Literal.Number.Bin */\n",
       ".mf { color: #666666 } /* Literal.Number.Float */\n",
       ".mh { color: #666666 } /* Literal.Number.Hex */\n",
       ".mi { color: #666666 } /* Literal.Number.Integer */\n",
       ".mo { color: #666666 } /* Literal.Number.Oct */\n",
       ".sb { color: #BA2121 } /* Literal.String.Backtick */\n",
       ".sc { color: #BA2121 } /* Literal.String.Char */\n",
       ".sd { color: #BA2121; font-style: italic } /* Literal.String.Doc */\n",
       ".s2 { color: #BA2121 } /* Literal.String.Double */\n",
       ".se { color: #BB6622; font-weight: bold } /* Literal.String.Escape */\n",
       ".sh { color: #BA2121 } /* Literal.String.Heredoc */\n",
       ".si { color: #BB6688; font-weight: bold } /* Literal.String.Interpol */\n",
       ".sx { color: #008000 } /* Literal.String.Other */\n",
       ".sr { color: #BB6688 } /* Literal.String.Regex */\n",
       ".s1 { color: #BA2121 } /* Literal.String.Single */\n",
       ".ss { color: #19177C } /* Literal.String.Symbol */\n",
       ".bp { color: #008000 } /* Name.Builtin.Pseudo */\n",
       ".vc { color: #19177C } /* Name.Variable.Class */\n",
       ".vg { color: #19177C } /* Name.Variable.Global */\n",
       ".vi { color: #19177C } /* Name.Variable.Instance */\n",
       ".il { color: #666666 } /* Literal.Number.Integer.Long */\n",
       "</style>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "importing Jupyter notebook from restaurants_data_cleaning.ipynb\n",
      "importing Jupyter notebook from review_data_process.ipynb\n",
      "importing Jupyter notebook from review_data_getData.ipynb\n",
      "importing Jupyter notebook from review_data_classify.ipynb\n",
      "importing Jupyter notebook from census2010_data_cleaning.ipynb\n",
      "importing Jupyter notebook from acs2013_data_cleaning.ipynb\n"
     ]
    }
   ],
   "source": [
    "import os.path\n",
    "import heapq\n",
    "from geopandas import GeoDataFrame, GeoSeries, read_file\n",
    "from geopandas.tools import sjoin\n",
    "from shapely.geometry import Point, Polygon\n",
    "import scipy.spatial as spatial\n",
    "from pandas import read_csv, to_numeric\n",
    "from math import log\n",
    "from numpy import where, concatenate, mean\n",
    "\n",
    "import jupynbimp\n",
    "import restaurants_data_cleaning, review_data_process, census2010_data_cleaning, acs2013_data_cleaning"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Source CRS: WGS 84/Lat,Long; Projected CRS: WGS 84/UTM Zone 12N\n",
    "CRS = {'GCS':'+init=epsg:4326', 'projected':'+init=epsg:32612'}\n",
    "USE_CACHE=True # if false, forces all modules to re-compute variables (takes significantly more time!) \n",
    "SQ_METERS_PER_SQ_KM=1000000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def normalizeByArea(geoDataFrame, column):\n",
    "    return geoDataFrame.apply(\n",
    "        lambda row: row[column]/row['area_sqkm'], axis=1\n",
    "    ) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Transform data into consolidated shapefile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def getData(fromCache=True):\n",
    "    # Retrieved 09-07-2016 from https://nz.yelp.com/dataset_challenge/dataset \n",
    "    dataDirectory = '../../data/'\n",
    "    outputName = 'restaurants_merge_clean.csv'\n",
    "    \n",
    "    if fromCache & os.path.isfile(dataDirectory + outputName ):\n",
    "        return read_csv(dataDirectory + outputName, header=0)\n",
    "    \n",
    "    else:\n",
    "        \n",
    "        # merge Yelp data\n",
    "        reviews = (review_data_process.getData(fromCache=fromCache)\n",
    "                   .set_index('business_id')\n",
    "                  )\n",
    "        \n",
    "        restaurants = (restaurants_data_cleaning.getData(fromCache=fromCache)\n",
    "                       .set_index('business_id')\n",
    "                       .join(reviews, how='inner', \n",
    "                             lsuffix='business_id', \n",
    "                             rsuffix='business_id')\n",
    "                       .rename(columns={'business_idbusiness_id':'business_id'})\n",
    "                      )\n",
    "        \n",
    "        # create spatial features\n",
    "        restaurants = GeoDataFrame(\n",
    "            restaurants,\n",
    "            geometry=[Point(x,y) \n",
    "                      for x, y in zip(restaurants['longitude'],\n",
    "                                      restaurants['latitude'])\n",
    "                     ],\n",
    "            crs=CRS['GCS']\n",
    "        )\n",
    "        \n",
    "        restaurants.to_crs(CRS['projected'], inplace=True)\n",
    "        restaurants.drop(['latitude', 'longitude'], axis=1, inplace=True)\n",
    "        \n",
    "        # Phoenix Convention Center in projected coordinates\n",
    "        phoenixCBD = GeoSeries(Point(400557, 3701492), crs=CRS['projected'])\n",
    "\n",
    "        restaurants['dist_CBD'] = restaurants.geometry.apply(\n",
    "            lambda point: phoenixCBD.distance(point)\n",
    "        )\n",
    "        \n",
    "        # Scottsdale shopping district in projected coordinates\n",
    "        scottsdale = GeoSeries(Point(413975, 3707095), crs=CRS['projected'])\n",
    "        \n",
    "        restaurants['dist_scottsdale'] = restaurants.geometry.apply(\n",
    "            lambda point: scottsdale.distance(point)\n",
    "        )\n",
    "        \n",
    "        # calculate distance from freeway exits\n",
    "        if fromCache & os.path.isfile(dataDirectory + 'shapefiles/motorway-exits/motorwayExits.shp'):\n",
    "            motorwayExits = read_file(dataDirectory + 'shapefiles/motorway-exits/motorwayExits.shp')\n",
    "            \n",
    "        else:\n",
    "            \n",
    "            # Retrieved 17/07/2016 from http://download.geofabrik.de/north-america/us/arizona.html\n",
    "            motorwayExits = read_file('/arizona-latest/roads.shp', \n",
    "                                     vfs='zip://' + dataDirectory + 'shapefiles/arizona-latest.zip'\n",
    "                                    ).to_crs(CRS['projected'])\n",
    "            \n",
    "            motorwayExits = (motorwayExits[motorwayExits['type'] == 'motorway_link'])\n",
    "            \n",
    "            # cache as shapefile\n",
    "            motorwayExits.to_file(dataDirectory + 'shapefiles/motorway-exits/motorwayExits.shp')\n",
    "            \n",
    "            \n",
    "        restaurants['dist_mwy_exit'] = restaurants.geometry.apply(\n",
    "           lambda point: motorwayExits.distance(point).min()\n",
    "        )\n",
    "        \n",
    "        \n",
    "        # for efficient nearest-neighbor queries\n",
    "        c2DTree = spatial.cKDTree(concatenate(restaurants.geometry.apply(lambda point: list(point.coords))))\n",
    "    \n",
    "        restaurants['nearest_neighbor_distance'] = restaurants.geometry.apply(\n",
    "           lambda point: c2DTree.query(list(point.coords), k=[2])[0][0][0]\n",
    "        )\n",
    "        \n",
    "        # competitor proximity = R statistic = r_obs / r_exp:\n",
    "        # r_obs = mean nearest neighbor distance of 5 nearest neighbors \n",
    "        # r_exp = 1/(2(n/study area)^0.5)(Lu & Wong, 2008)\n",
    "        r_exp = 1 / (2 * (restaurants.shape[0] / restaurants.geometry.unary_union.envelope.area)**0.5)\n",
    "        restaurants['competitor_proximity'] = restaurants.geometry.apply(\n",
    "           lambda point: (restaurants.iloc\n",
    "                          [c2DTree.query(list(point.coords), k=[2, 3, 4, 5, 6])[1][0]]\n",
    "                          ['nearest_neighbor_distance']\n",
    "                          .mean() / r_exp\n",
    "                         )\n",
    "        ) \n",
    "        \n",
    "        # merge demographic (US census) data\n",
    "        census2010 = census2010_data_cleaning.getData(fromCache=USE_CACHE)\n",
    "\n",
    "        # Retrieved 10-07-2016 from https://www.census.gov/cgi-bin/geo/shapefiles/index.php?year=2010&layergroup=Block+Groups\n",
    "        blockGroups2010 = read_file('/tl_2010_04_bg10.shp', \n",
    "                                    vfs='zip://' + dataDirectory + 'shapefiles/tl_2010_04_bg10.zip')\n",
    "        blockGroups2010.GEOID10 = to_numeric(blockGroups2010.GEOID10)\n",
    "        blockGroups2010 = (blockGroups2010[['GEOID10','geometry']]\n",
    "                           .rename(columns={'GEOID10':'GEOID'})\n",
    "                           .to_crs(CRS['projected'])\n",
    "                          )\n",
    "        blockGroups2010['area_sqkm'] = blockGroups2010.geometry.area/SQ_METERS_PER_SQ_KM \n",
    "        blockGroups2010 = blockGroups2010.merge(census2010, on='GEOID')\n",
    "\n",
    "        # Normalize count data by block group area\n",
    "        blockGroups2010['population_density'] = normalizeByArea(blockGroups2010,'population_total')\n",
    "        blockGroups2010['home_mortgage_density'] = normalizeByArea(blockGroups2010,'home_mortgages')\n",
    "        blockGroups2010['home_owner_density'] = normalizeByArea(blockGroups2010,'home_owners')\n",
    "        blockGroups2010['renter_density'] = normalizeByArea(blockGroups2010,'renters')\n",
    "        blockGroups2010['household_density'] = normalizeByArea(blockGroups2010,'total_households')\n",
    "        blockGroups2010['family_household_density'] = normalizeByArea(blockGroups2010,'family_households')\n",
    "        blockGroups2010['single_household_density'] = normalizeByArea(blockGroups2010,'single_households')\n",
    "        blockGroups2010['hispanic_latino_population_density'] = normalizeByArea(blockGroups2010,'population_hispanic_latino')\n",
    "        blockGroups2010['white_population_density'] = normalizeByArea(blockGroups2010,'population_white')\n",
    "        blockGroups2010['black_population_density'] = normalizeByArea(blockGroups2010,'population_black')\n",
    "        blockGroups2010['native_american_population_density'] = normalizeByArea(blockGroups2010,'population_native_american')\n",
    "        blockGroups2010['asian_population_density'] = normalizeByArea(blockGroups2010,'population_asian')\n",
    "\n",
    "        blockGroups2010.drop(['population_total',\n",
    "                              'home_mortgages',\n",
    "                              'home_owners',\n",
    "                              'renters',\n",
    "                              'total_households',\n",
    "                              'family_households',\n",
    "                              'single_households','population_hispanic_latino',\n",
    "                              'population_white',\n",
    "                              'population_black',\n",
    "                              'population_native_american',\n",
    "                              'population_asian'\n",
    "                             ], \n",
    "                             axis=1, inplace=True\n",
    "                            )\n",
    "        \n",
    "        acs2013 = acs2013_data_cleaning.getData(fromCache=False) \n",
    "\n",
    "        # Retrieved 10-07-2016 https://www.census.gov/cgi-bin/geo/shapefiles/index.php?year=2013&layergroup=Block+Groups\n",
    "        blockGroups2013 = read_file('/tl_2013_04_bg.shp', \n",
    "                                    vfs='zip://' + dataDirectory + 'shapefiles/tl_2013_04_bg.zip')\n",
    "        blockGroups2013.GEOID = to_numeric(blockGroups2010.GEOID)\n",
    "        blockGroups2013 = (blockGroups2013[['GEOID','geometry']]\n",
    "                           .to_crs(CRS['projected'])\n",
    "                          )\n",
    "\n",
    "        blockGroups2013['area_sqkm'] = blockGroups2013.geometry.area/SQ_METERS_PER_SQ_KM \n",
    "        blockGroups2013 = blockGroups2013.merge(acs2013, on='GEOID')\n",
    "\n",
    "        # Normalize count data by block group area\n",
    "        blockGroups2013['density_education_highschool'] = normalizeByArea(blockGroups2013,'education_highschool')\n",
    "        blockGroups2013['density_education_undergraduate'] = normalizeByArea(blockGroups2013,'education_undergraduate')\n",
    "        blockGroups2013['density_education_postgraduate'] = normalizeByArea(blockGroups2013,'education_postgraduate')\n",
    "\n",
    "        blockGroups2013.drop(['education_highschool',\n",
    "                              'education_undergraduate',\n",
    "                              'education_postgraduate'\n",
    "                             ], \n",
    "                             axis=1, inplace=True\n",
    "                            )\n",
    "        \n",
    "        # spatial join restaurants with census block groups\n",
    "        restaurants.reset_index(inplace=True)\n",
    "        restaurants_census2010 = (sjoin(restaurants[['business_id', 'geometry']], blockGroups2010, how='inner')\n",
    "                                  .drop(['geometry', 'index_right', 'GEOID', 'area_sqkm'], axis=1)\n",
    "                                 )\n",
    "        restaurants_acs2013 = (sjoin(restaurants[['business_id', 'geometry']], blockGroups2013, how='inner')\n",
    "                               .drop(['geometry', 'index_right', 'GEOID', 'area_sqkm'], axis=1)\n",
    "                              )\n",
    "        restaurants = GeoDataFrame(restaurants_census2010.merge(restaurants_acs2013)\n",
    "                                   .merge(restaurants), \n",
    "                                   crs=CRS['projected']\n",
    "                                  )\n",
    "  \n",
    "        return restaurants"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
